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Abstract 



The perfect quenching of spin tunneling first predicted for a model with 
biaxial symmetry, and recently observed in the magnetic molecule Fes j is fur- 
ther studied using the discrete phase integral (or Wentzel-Kramers-Brillouin) 
method. The analysis of the previous paper is extended to the case where the 
magnetic field has both hard and easy components, so that the Hamiltonian 
has no obvious symmetry. Herring's formula is now inapplicable, so the prob- 
lem is solved by finding the wavefunction and using connection formulas at 
every turning point. A general formula for the energy surface in the vicinity 
of the diabolo is obtained in this way. This formula gives the tunneling apm- 
plitude between two wells unrelated by symmetry in terms of a small number 
of action integrals, and appears to be generally valid, even for problems where 
the recursion contains more than five terms. Explicit results are obtained for 
the diabolical points in the model for Fes- These results exactly parallel the 
experimental observations. It is found that the leading semiclassical results 
for the diabolical points appear to be exact, and the points themselves lie on 
a perfect centered rectangular lattice in the magnetic field space. A variety 
of evidence in favor of this perfect lattice hypothesis is presented. 

75.10Dg, 03.65.Sq, 75.50Xx, 75.45.+j 
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I. INTRODUCTION 



A. The story so far 

In a previous paper |J, hereafter cited as I, we studied the tunneling of a spin governed 
by the Hamiltonian 

H = -k 2 J 2 z + (h - k 2 )Jl - gf i B J ■ H, (1.1) 

where J is a dimensionless spin operator, and k\ > k 2 > 0. This Hamiltonian is the simplest 
descriptor of the magnetic properties of the molecule [(tacn) 6 Fe802(OH)i 2 ] 8+ (or just Feg 
for short), with J = 10, k\ w 0.33 K, and k 2 « 0.22 K [0-f|]. Interest in this molecule arises 
because of its rich low temperature magnetic behavior, which include hysteresis at the level 
of one molecule |J, and more recently, the discovery || of an entire lattice of diabolical 
points in its magnetic spectrum, a subset of which was predicted to exist earlier []8| . It 
is the latter property that we wish to continue investigating in this paper. 

In paper I, attention was confined to the case where the external magnetic field H is along 
the hard axis x. In this case, the classical energy, which may be viewed as the expectation 
value (TC) of the Hamiltonian in a spin-coherent state, is symmetric about the xy plane, 
and the problem is analogous to that of a massive particle in one dimension in a reflection 
symmetric double well potential. It is then natural to consider the tunneling between the 
symmetrically related states localized in the left- and right-hand wells. In the magnetic case, 
the analogous states are those with predominatly positive and negative values of J z , at least 
as long as H x is not so large as to bring the classical minima very close to the x axis. The 
surprise is that the tunnel splitting between the ground states oscillates as a function of H x , 
vanishing exactly at a series of points. In fact, the splitting between higher pairs of levels 
also vanishes at just these points, as noted in I and earlier studies P,|10 . 



In addition to observing the quenching of ground state tunneling when H||x, however, 
Wernsdorfer and Sessoli || also performed experiments with H z ^ 0. The reflection sym- 
metry of the classical energy is now lost, but if the value of H z is chosen properly, it is 
possible to bring an excited state in the positive J z well into approximate degeneracy with 
the ground state of the negative J z well. The new discovery by them is that if H x is now 
varied, the tunnel splitting between the degenerate or quasi-degenerate levels again oscil- 
lates. It is theoretically understood that if both H z and H x are properly tuned, the splitting 
vanishes exactly in this case too P,|TT] . Experimentally, of course, one can never see a perfect 



zero in the splitting, and Wernsdorfer and Sessoli only see a minimum in the Landau-Zener- 
Stiickelberg transition rate between the levels in question. The minima are so deep, however, 
that there is little doubt that the underlying tunneling matrix element is quenched. 

When H||x, the Hamiltonian (|1.1| ) is invariant under a 180° rotation about x, and the 
quenchings can be understood from the viewpoint of the von Neumann- Wigner theorem 
as allowed crossings of energy levels with different parities under this rotation. A similar 
argument can be made when H||z. When H has both x and z components, however, the 
Hamiltonian has no obvious symmetry, and the above theorem states that an intersection 
of two energy levels is infinitely unlikely as a single parameter in the Hamiltonian is varied. 



For a real symmetric Hamiltonian, it is known [12| that one must vary two parameters to 



obtain an intersection. Since we can choose the matrices of both J, and J T to be real in the 
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J z basis, these two parameters can be taken as H x and H z . The isolated points ]T3[ in the 
H x -H z plane where any two levels intersect are precisely what Herzberg and Longuet-Higgins 
H call conical intersections and what Berry and Wilkinson [7fl call diabolical points. The 
latter terminology originates in the resemblance of the energy surface — a double cone with 
a common vertex — to an Italian toy called the diavolo. 



B. Content and plan of this paper; the perfect lattice hypothesis 

In this paper we shall allow H to lie in the x — z plane, with a view to studying the 
tunneling in the asymmetrical well, and locating the diabolical points. As in paper I, our 
analysis is based on the discrete phase integral (DPI), or Wentzel-Kramers-Brillouin method 
|15| . This method is semiclassical in character, with 1/ J playing the same role as h in the 



continuum phase integral method. To introduce this method, let us write the Schrodinger 
equation 7i\ip) = E\ip) in the J z eigenbasis With J z \m) = m\m), (m\ip) = C m , 

(m\TC\m) = w m , and (m\H\m') = t m>m i (m ^ m'), we have 

m+2 

^ ] tm,n,C- n ~\~ W m C m EC mi (1-2) 
n=ra— 2 

where the prime on the sum indicates omission of the n = m term. A vivid picture of the 
approximation can be obtained if we think of Eq. (|1.2| ) as the tight-binding model for an 
electron in a one-dimensional lattice with sites labelled by m, and on-site (w m ) and hopping 
(t m m±1 , t m ^ m ±2) energies. If J 3> 1, these matrix elements vary slowly with m, on a length 
scale of order J in fact. We may then use the approximation of semiclassical dynamics 
by working entirely in terms of wavepackets whose spatial extent is much less than the 
length scale over which the properties of the lattice vary, i.e., J, and whose spread in Bloch 
vectors is much less than the width of the Brillouin zone, i.e., 2tt. These ideas have close 
counterparts in the continuum quasiclassical method, and the DPI method is nothing but 
the discrete analog. 

When H||x, the problem can also be approached using instantons — indeed the oscillations 
in the splitting were discovered in this way. When H has other components besides H x , 
however, the DPI method is, to our knowledge, the only successful one to date. Villain and 
Fort's approach M is also an approximate version of this method that makes artful use of 
some special features of the Fe 8 problem and works for small values of the field. Our analysis 
is more prosaic, and almost self-evident once one has understood how to deal with the new 
feature in Eq. ( |1.2| ) — the presence of second neighbour hopping. This gives rise to novel 
turning points with no continuum analogues. Connection formulas for these turning points 



are given in Ref. [ IB[ . We have quoted the results of our analysis before [O,0(b), but the 
details are only being presented here. 

Our main result for the specific Hamiltonian ( |1 . 1[ ) is for the locations of the diabolical 
points. We find that the £'th level in the negative J z well (where £' = denotes the lowest 
level) and the £"th level in the positive J z well are degenerate when (see Fig. [I]) 

H z (£',£") _ VX(£" - g) 
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^n = yHlA [j- n -Ui'+i" + l)}, (1-4) 

with n = 0, 1, . . . , 2J - (£' + I" + 1). Here, A = k 2 /k 1 , and F c = 2k 1 J/g f i B . 

It should be stressed that Eqs. (|1.3|) and (|1.4| ) only represent the first terms of a series in 
1/J, and that our DPI calculations give no reason to believe that the higher order terms are 
absent. Yet a large amount of empirical evidence suggests just this, i.e., that the formulas 
are exact as written] That the diabolical points lie on a perfect centered rectangular lattice, 
and that many pairs of levels are simultaneously degenerate, we refer to as the perfect lattice 
hypothesis. It has previously been made for the ground state on the H z — line by Villain 
and Fort 0, and extended to include simultaneous degeneracy of the higher states by us 
|l(|(b). This was shown by perturbation theory in A to 0(A 3 ) for all J, and analytically for 
J < 2. We have since found analytically that for J = 3/2, the additional diabolical point 
at (H X ,H Z ) = (VT^A, VA)# C / 3 is exact - 0f 

course, Kramers theorem guarantees double 
degeneracy of all energy eigenvalues at the point H x — H z — for all half-integral J. We 
have also checked this hypothesis numerically for a variety of values of A, and J up to 10. 
We present our results for J = 5 in Table I. The deviations in the locations of the diabolical 
points from the perfect lattice hypothesis predictions are never more than 10~ 10 , certainly 
much less than order 1/J. In fact, the values listed are below the numerical tolerance that 
we prescribed. 

Further support for the perfect lattice hypothesis comes from the fact that it is consistent 
with the following duality property of the Hamiltonian ( |1 . 1| ) . If we set k\ — 1, Ti. may be 
written as H(X, H x , H z ), showing its dependence on the three parameters, A, H x and H z . 
By a 90° rotation about the axis (x + z)/ v2, we obtain the transformation 

H(\, H x , H z ) <- -H{\ - A, H z , H x ). (1.5) 

In particular, the spectra of the two Hamiltonians are so related, and ranking the levels is 
order of increasing energy, we see that if the levels with ordinal numbers k and k + 1 are 
degenerate when H x = f x (X) and H y = f y (\), then level numbers 2 J + 2 — k and 2 J + 1 — k 
are degenerate when H x = f y (l — A) and H y = f x (l — A). These conditions do not constrain 
the functions f x and f y in any real way, however, so the discovery of formulas (|1.3|) and ( |1.4|) 
must presently be put down to serendipity. 

The plan of our paper is as follows. In Sec. [0], we analyze the asymmetric double well 
problem in completely general terms. We briefly review the DPI method and the asymmetric 
double-well problem in one dimensional quantum mechanics. In contrast to the situation 
that prevails in that case, we must deal with four independent DPI wavefunctions at every 
point. Consequently the matching problem is harder and more subtle. Its study forms the 
bulk (subsections C to G) of the analysis. The main result of this analysis is a formula for 
the diabolo [Eq. (|2.63|) ], and a two level system Hamiltonian [Eq. (|2.64j )l. The tunneling 
amplitude appearing in these formulas is found in a general form [Eq. (|2.60|) 1 involving a 
small number of action integrals running between various turning points. Readers who are 
not interested in the details of the analysis should skip to subsection H and read on from 
Eq. ( |2.63|) . They may also find the last paragraph of Sec. [H] interesting. 

In Sec. PH we apply the general analysis to the model ( |1 . 1|) for Fes- To keep the analysis 
tractable, we will assume that H z /H c is small, although there is no reason why one could 
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not apply the results of Sec. |T| numerically for arbitrary values of H z . It may in fact 
be interesting to do this to further investigate the perfect lattice hypothesis. Further, as 
discussed in Sec. [II], the duality property of TC immediately extends our work to large H z 
and small H T . 



II. GENERAL FORMULA FOR DIABOLO IN TERMS OF ACTION INTEGRALS 

A. Summary of DPI method, critical curves, etc. 

Very briefly, the DPI analysis proceeds as follows. (See Ref. |16| and paper I for details.) 
Let the energy of an electron wavepacket in the equivalent tight-binding model be given 
by 7i sc (q,m) where q and m are the mean wavevector and position of the wavepacket, 
respectively. Holding m fixed, we may think of 7i sc {q, m) as a dispersion relation for the local 
energy band. The mean velocity of the wavepacket is then given by v (q, m) = &H SC (q, m) / dq. 
If we regard m as a continuous variable, and approximate w m , t m , TO ±i, and t m , m ±2 by smooth 
functions w(m), ti(m), and t2(m), then 

TC sc (q, m) = w{m) + 2t 1 (m) cosg + 2t 2 (m) cos(2g), (2.1) 
v(q,m) = — 2sing(m)(ti(m) + 4t 2 (m) cosg(m)). (2.2) 

In close analogy with the continuum method, one can show that the quasiclassical solution 
to the wavefunction for a given energy E is given by linear combinations of solutions of the 
form 

C m ~ exp (i f q(m')dm'^\ , (2.3) 



>v{m) 

where q(m) and v{m) are given by 



E = H sc (q(m) , m) , v (m) = v (q(m), m). (2.4) 



It follows from Eq. ( p.3|) that the quasiclassical solution breaks down when v{m) = 0— 
which we call a turning point, as in the continuum case. At these points, the solutions must 
be augmented by connection formulas. In addition to the values q = and q = ir, v{m) also 
vanishes at q — q*(m) = cos~ 1 (—ti(m)/4t2(m)). We will get turning points whenever E 
equals Uo(m), U n (m), or U*(m), where these three functions are ?i sc (0, m), H sc (vr,m), and 
7i sc (q*,m) respectively. These three curves, which we call critical curves, collectively play 
the same role as the potential energy in the continuum case. 

For Fe 8 , in the same reduced variables as in I = m/J, energies in units of kiJ 2 , 
J = ( J + |)], the on-site and hopping energies are given by 

w(m) = i(l + A)(l-/i 2 ) -2h g fi, (2.5) 
t 1 (m) = -h x (l- f i 2 ) 1 ^, (2.6) 
t 2 (m) = I(l-A)(l-/i 2 ). (2.7) 

Here, A = k 2 /ki, h x = JH X /JH C , h z = JH Z /JH C , and H c = 2k 1 J/gfi B . Thus the critical 
curves are given by 
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fi 2 - 2h x (l - \l 



2h x (l 



2)1/2 
2x1/2 



2h z fi, 
2h z fi, 



UM = A(l 



M 2 )- 



h 2 



(2.8) 
(2.9) 

(2.10) 



These curves are shown in Fig. |2|. The energy U^{m) is the upper edge of the band TC sc (q, m) 
for all m. The lower band edge is given by U (m) for values of m close to ±J, and by U*(m) 
for m in the central region. In the central region, the second neighbour hopping element 
t2(m) is sufficiently large that the local energy band H, sc (q,m) has its global minimum not 
at q — 0, but at q*. Since the energy U lies in the band, U n > U > U*. By contrast, in the 
outer m regions, Ti. sc (q,m) has only one minimum (at q = 0) and only one maximum (at 
q = 7r) for real g. The energy £7* lies outside the band, so that once again > U > U*. 
The curves Uo(m) and U*(m) touch one another with a common tangent at m — ±m*, 
where 



J 



h 2 . 



1/2 



(2.11) 



The turning point where E — U* is special if it happens at a value of m where U* lies 
below the lower band edge. The wavevector q* at the turning point is then complex, and 
the wavefunction C m changes from a decaying (or growing) exponential on one side to a 
decaying (or growing) exponential with an oscillatory envelope on the other side. These 
turning points are the new feature caused by second neightbor hopping that we referred to 
in Sec. |. 



B. Nature of asymmetric double well wavefunctions 

To understand how the DPI solutions are to be used to find the eigenstates, it is useful 
to discuss the corresponding problem for an asymmetric double well in the continuum case. 
Suppose the potential is as drawn in Fig. [3]. The potential minima are at Xq±, and for 
the energy E drawn, the turning points are at x' a and x' b in the left well, and x'l and x" a 
in the right well. We use the quasiclassical approximation to find a wavefunction ip'(x) on 
the left hand side as follows. First, we choose the solution which decays exponentially as 
x — x' a — ► — oo. This solution is matched via connection formulas at x' a to an oscillatory 
solution in the region x' a < x < x' b . We then use connection formulas at x' b to find the 
quasiclassical solution in the region x > x' b . We repeat this procedure on the right hand 
side to find a wavefunction ip"(x) that decays to zero as x — x" —>■ oo. The last step is to 
demand that the wavefunctions ip'(x) and ip"(x) be the same in the central region, i.e., in 
the vicinity of x = 0. This demand will be unsatisfiable for an arbitrarily chosen energy E, 
and will provide one with the eignevalue condition. 

There are two remarks that we wish make about the above procedure. The first remark 
concerns the basic nature of the solution. In general, in the central region, the left solution 
ip'{x) will be a linear combination of two parts, ip' d {x) and ip'gix), that are exponentially 
decaying and growing as x — x' b increases, respectively. Likewise the right solution, tp", will 
be a sum of parts that decay (ip'^) and grow (ipg) as x b — x increases. The key point is that 
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the growing parts ip' and ip" must be present in an eigenstate, for without them, there is no 
way that the values and slopes of if)' and if)" could be made to agree at x = 0, say. 

The second remark is technical. If the potential well is reasonably parabolic near the 
minimum at xq-, the Schrodinger equation can be solved directly for any choice of E in 
terms of parabolic cylinder functions, and we can always find a linear combination that will 
decay to zero as x — x' a — > — oo. This linear combination will have both growing and decaying 
pieces as x — x' b grows. In this way we can obtain the wavefucnction ip'{x) on the entire 
left hand side without using connection formulas at x' a or x' b . Once one has found ip'(x) 
sufficiently far to the right of x' b in this way, one can write the parabolic cylinder functions 
in quasiclassical form which can then be extended in this form all the way to x — 0. The 
right-hand wavefunction can be treated in the same way. This device leads to considerable 
savings in labor. 

We now apply these ideas to our problem. In what follows, we will denote quantities 
pertaining to the left hand solution or the left hand side of the well by either a single prime 
or a suffix — , while analogous right-hand quantities will carry a double prime or a + suffix. 
We consider an energy E as drawn in Fig. |], which leads to turning points at m' a , m' b , m' c 
on the left hand side, and m", m b , and vn!' a on the right. The wavefunction C' m will decay 
away from the well bottom as m — m' a decreases, oscillate in the classically allowed region 
m' a < m < m' b |fT7| , [IS[| . In the region just to the right of m' b it will consist of a decaying 
part and a growing part. The new feature will be encountered at m' c where E = U*. For 
m > m' c , both the growing and decaying parts will acquire oscillatory envelopes. Similar 
remarks apply to the right side wavefunction C'^. 



C. DPI wavefunction in the leftmost forbidden region 

We are now ready to find the wavefunction explicitly. Let us start constructing C' m from 
the left. The Hamilton- Jacobi equation in Eq. (|2.4j) has the general solution 



, , -him) ± \t 2 (m) - U 2 (m)f(m)] 1 / 2 

cosgm = u ; ^Jl^Jl . 2.12 

U 2 (m) 

f(m) = w(m) — 2t 2 (m) — E. (2.13) 

This leads to four values of q(m) for any E, since if q is a solution, so is —q. In the region 
m ^ m ' a i since E < Uq, all four solutions are pure imaginary. We write the two which lead 
to decaying wavefunctions as \i — n' a becomes large and negative as 

9i )2 = -^i, 2 (/i), (2.14) 
where k 2 > Ki > 0, and the corresponding DPI solutions as 

C' m l = A'\v 1 {m)\- 1 ' 2 e-ycg U J™ qi(m')dm'] , (2.15) 

C' m ,2 = B'\v 2 (m)\- 1/2 exp U J™ q 2 (m')dm^j , (2.16) 

with Vi(m) = v(qi(m),m). We take A' and B' to be real without any loss of generality. Note 
that we have left the lower limit of the phase integral for C' m l unspecified and written it as 
m' c for C' m2 . The reasons for this will become clear shortly. 
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To see how these two solutions behave as m approaches m' a , and continues beyond this 
point, let us note that 



cosh k x 2 = 1 1 ^ 1 1 t — . 2.17 

As m — > m' a —, cosh/ti — > 1, i.e., Ki — > 0, while coshfi^ — > — 1 + 2|ti|/4t2 > 1- As we cross 
the point m' a , q\ will become real, while g 2 will continue to be pure imaginary and large. 



Thus the solution C' m2 continues to hold at m' a , while C' ml breaks down at m' a fll8| , and 



must be related to a the solution for m > m' a by a connection formula. It is clear that the 
wavevector(s) for C' m2 will continue to be given by Eq. ( |2.17| ) as m — > m' b —, while those 
for C' m 1 will again approach zero, necessitating the use of connection formulas to go on 
to m > m' b . It is here that the technical remark about sidestepping the use of connection 
formulas that was made in connection with the continuum antisymmetric double-well is 
relevant. The solution C' ml can be approximated by a harmonic oscillator wavefunction 
provided the energy E is not very far from the minima of Uq. The asymptotic forms of this 
wavefunction give us the quasiclassical wavefunction in the regions m < m' a and m' b more 
simply. We therefore turn to this subproblem. 

D. Jumping across the potential well 

The assumption that E — Uq(hq±) is not very large, means that qi(m) is never very far 
from zero, and we may expand Ti. sc (q, m) in powers of q and m + m' Q . As in I, we write 

H sc ( 9) m) = £_ + ^- 9 2 + ^M„w!(m - m -,) 2 + • • • ■ (2.18) 

Since q and m are conjugate variables, we can write C' m l as the solution to the differential 
equation 

H sc (-idm,m)C' mjl = EC' m!l . (2.19) 

Introducing two new variables z and v' by the equations 

m = m Q _ + (2M_uo _y 1/2 z, (2.20) 
E = E„ + (2.21) 



within the approximation (|2.18|) , the differential equation becomes that for the parabolic 
cylinder functions: 



d 2 ( . 1 z 2 
dz- 2 + V + 2"T 



C' = 0. (2.22) 



If we take as the two linearly independent solutions the standard forms D v i(z) and D v i{—z) 
fHfl , the former must be rejected as it diverges for z — > — oo. We accordingly write 

C,i = ^'(-l)%(-4 (2-23) 
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where A' is the constant in Eq. (|2.15| ), and the additional factor (— l) e ' , where we shall define 
£' shortly, is another constant introduced for later convenience. 

As z — > — oo, D v /(—z) ~ (— z) u ' e~ z / 4 , and one can show with a little work that [modulo 
the factor (— 1)^] this is indeed the DPI form for C' ml with the approximation ( 2.18 ) [20 
For z —>■ +oo, on the other hand, 



DA-z) ~ cos(7rz/')^'e- 22 / 4 1 



v\v' 



2z 2 



+ ^-v-( 1 + <^^ + ...). ( , 24) 

Note that this form has both decaying and growing components. In fact, the latter compo- 
nent vanishes only if v' is a positive a integer. As explained earlier, it is essential for our 
DPI solution C' m to contain a growing component. We therefore allow for its presence by 
writing 

v' = £' + —, (2.25) 
No- 
where IS is a positive integer, and e' is a shift defined to lie in the interval (—1/2, l/2)o;o_. 
In fact, we expect e' to be very close to zero for any state in which there is large probaility 
of finding the particle in the left well. We can make the vanishing of the growing component 
in C' m 1 more manifest by writing 

1 Siu( ""\(l + i/) « (-if'+V'O — • (2-26) 



T(-U') 71 LO 

Combining Eqs. (|2.23|) -( 2.26 ), we thus find that for m beyond m' b , 



C' 1 « A'\ cos — z v e~ z /4 - V27t(£' !)— z' v - 1 e 2 /4 . (2.27) 



E. DPI form in ordinary forbidden region 



The next step is to write the solution for C' m in such a way that it holds in the entire 



region m' b < m < m' c [Pjfl . For (7^2, Eq. ( p.!6[ ) already meets this demand, since just as at 
the turning point m' a , q 2 stays imaginary and negative as m passes through m' c . For C' m 1( on 
the other hand, Eq. fl2.27| ) only holds in a region where the parabolic approximation to Uq 
is good, and may not hold near m' c . It is clear, however, that the wavevector qi associated 



with C' n 
write CL 



m,l 



m,l 



is again given by —%K\ with K\ given by Eq. 
in the DPI form 



( |2-17| ). Hence it must be possible to 



CL 



m,l 



\vi{m)\ 



-1/2 



Q' exp I — i j qi (m')dm'\ + R' exp lij q 1 (m')dm' 



(2.28) 
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where Q' and R' are coefficients which we expect to be proportional to A'. To find these, 

let us first calculate the phase integral in the parabolic approximation. Since qi = at m' b , 
we have 

E = H sc (0,m' b ) =H sc (tK l (m) } m), (m > m' b ). (2.29) 
Using Eq. fl2.18| ), we get 

/ti(m) « M-W Q -[{m - m -) 2 - (m' 6 - m _) 2 ) 1/2 ]. (2.30) 
Therefore, if (m — m' b ) 3> (mj, — mo-), we obtain 



l 



V / \ m b~ m o- ) 

x exp f--M_u; _(m - m _) 2 ) , (2.31) 

where we have used the fact that 

\M_ul_(m' b - m _) 2 = W + i) Wo _. (2.32) 



Next we note that (a) from the definition of z, Eq. (|2.20|) , the last exponential in Eq. (|2.31|) 



is nothing but exp(— z 2 /4), and that (b) |t>i(m)| « /ti (m)/M_ c<Jo-(m — mo-). Since 
m — m _ oc it follows that the firat term in Eq. ( [2.28|) varies as z v exp(— z 2 /4), i.e., 
precisely as the first term in Eq. ( ^.27| ). Thus Q' is indeed proportional to A', and a little 
algebra plus the use of Eq. ( |2.32|) gives 

Q'=(^-] V4 e-W + I) W + \)W + D cos ^A> = a' A'. (2.33) 
In the same way, we can show that 



1/4 



2M..J V 2; 1 W too- 



The definitions of the factors a' and /?' which we have introduced for later convenience can 
be read off these equations. Note that since we took A' to be real, Q' and R' are also real. 

To summarize where we are, the complete DPI solution for C' m in the region m! b < m < m' c 
is given by the sum of Eqs. ( |2.16|) and ( |2.28| ). The terms in B' and R' are exponentially 
growing with increasing m, while the term in Q' is exponentially decreasing. The next step 
is to connect this solution to the DPI form in the region m' c < m. As already stated, the 
turning point m' c is the irregular one under the barrier, where exponentially growing and 
decaying wavefunctions acquire oscillatory envelopes as it is crossed. 

Before we use the connection formulas at m' c , it is useful to see how the quasiclassical 
wavevector behaves near this point. Since cosg(m^,) = cosg* = — ti(m^)/4t 2 (m[,), it follows 



that the discriminant in Eq. (|2.12|) vanishes at m' c , and both qi and qi tend to the same value 
— in' , where 
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cosW c = (|ti|/4i 2 ) m=m ' c . (2.35) 

As we cross m' c , the discriminant in Eq. ( |2.12| ) becomes negative and cosg(m) [and therefore 
q{m)} becomes complex. We separate q(m) into its real and imaginary parts, and write two 
distinct solutions as 

q^gijn) = ±i«(m) + ( m > m ' c ) (2.36) 

where both k and \ are real and positive. The subscripts d and g stand for 'decaying' and 
'growing'. 



F. DPI form in oscillatory forbidden region 



The connection formulas to be used at m' c were derived in Sec. IV of Ref. [ 16| . The 
parts of C' m multiplying B', Q', and R! correspond respectively to the cases there labelled 
(o"i, a 2) = (+1, —1), (—1, +1), and (—1, —1). The B' part is given by 



C' m ,2 = B' 



Sg(m) 



7T 



exp i / q g {m')dm' + 



c.c. 



where 



s g {m) = 8t2{m) sinh K(m) sinx(m) sinq g (m) 

= 8^2 sinh k sin x(si n X cosh k — i cos x sinh k). 

Likewise, the coefficient of Q', which we call part la, connects to 

1 



Cm, la — e Q 



sAm) 



exp % \ q d (m')dm' + c.c. 



where 



Sd{m) = 8t2(m) sinh k(iti) sin x(m) singd(m) = [s 9 (m)]*, 
and T' is the phase integral, 



K(m')dm', 



(2.37) 



(2.38) 



(2.39) 



(2.40) 



(2.41) 



which we acquire in changing the lower limits of the m integrals from m' b to m' c . Lastly, the 
term in R', which we call part lb, connects to 



C - -p f ' Ft' 



exp (i I q g {m')dm'\ + c.c. 
a (m) ' V Jm 'c J 



[2A2) 



Equations ( p.37| ), (|2.39|) , and ( |2.42|) , give us the complete wavefunction C' m in the central 
region near m = 0. To simplify the writing, we denote 
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[AiK(m') + iA2x( m/ )]^ m ') 



(2.43) 



where Ai and A2 can be ±1 independently. In other words, the subscripts on $' give the 
signs of the real and imaginary parts [$' ++ is the integral of (n + ix), $'_ + that of (— K + ix), 
etc.]. The complete DPI solution for C' m can then be written as 



Q 



\Js* g (m) 



+ (\e T 'R' - iB' 



(m) 



\fsgim) 



+ c.c. 



(2.44) 



The wavefunction from the right, can now be written down at once. We define 
quantities with double primes in exact correspondence with those for C' m . The analog of 
Eq. ( [Q4D is then 

/-in 



e- r "Q" 



s*(m) 



+ (±e r "i?" - z£?" 



e *++M 



+ c.c. 



(2.45) 



The only issue requiring any thought is what the sign suffixes in $" should mean. By defining 
a new variable n = — m, so that the problem for becomes completely isomorphic to that 
for C' and then transforming back to m, one can show that 



r" 



$" + (m) 



K(m')dm' , 

[— Ac(m') + ix{m')]dm' , 



(2.46) 
(2.47) 



etc. Note that since m < m" in the center, and m" < m'l, these integrals are written so that 
the lower limit is less than the upper limit. Thus the suffixes on give the true signs of 
its real and imaginary parts. 



G. Matching of left and right wavefunctions 

It remains to see if Eqs. ( |2.44j ) and ( |2.45| ) are the same function. We note that $" + (m) 

has the same integrand as $^ (m) if m is taken to be the upper limit for both integrals. We 

further note that the remaining m dependence in both terms is (s*) -1 / 2 . Similar remarks 
apply to the and <3>'_j__ terms. Thus, we conclude that C' m will equal if the following 
conditions are obeyed: 



e 



-r"Q" e *'l + M = (l e v R > + iB >) e *;-M (2.48) 
'Q' e «'_+M = (l e T "R" + iB") e *+-H. (2.49) 
To simplify these conditions we note that 



$' + _(m) — $" + (m) = / [n(m) — ix{m)} dm, 

= T c -tA c , (2.50) 



12 



where T c and A c are the real and imaginary parts of the integral, and the subscript 'c' 
indicates that the integrals extend over the central region of m. Equations (gjgp and (EOSD 
can now be written as 



r c -iA c 



Q"= (l e r 'R f + iB f ) e r "e 
Q'= he T "R" + iB») e r 'e r - lA 



(2.51) 
(2.52) 

If we recall that [see Eqs. ( |2.33| ) and ( |2.34| )1 Q' and R! are proportional to A', and likewise 
for Q" and R", these equations are two complex equations in the four real quantities A', B', 
A" and B". To solve them we first note that the imaginary parts on the right hand sides 
must vanish. This yields 



B' 
B" 



\R'e r 'tan A c , 
\R"e T " tanA c . 



Substituting these in Eqs. ( |2.51|) and ( |2.52| ), we obtain after some simplification 

R' = 2e- Ta cosA c Q", 
R" = 2e- v ° cosA c Q', 

where To is the total Gamow factor 



G 



Ki(m)dm, 



(2.53) 
(2.54) 



(2.55) 
(2.56) 



(2.57) 



and the subscript '1' on k is to remind us that we must use the imaginary part of the 
wavevector that goes to zero at the turning points m' b and m'l. 



H. The eigenvalue condition, and the diabolo 



The simplest way of solving Eqs. (|2.55|) and ( |2.56| ) is terms of the ratios a' and f3' defined 



in Eqs. ( |2.33| ) and ( |2.34j ), and the analogous ratios a" and (3". Equating the products of the 
left hand and right hand sides, and simplifying a little, we get 

e ' e " = 4 ^7^o-^o+ cos 2 A c e~ 2TG . (2.58) 

This is our eigenvalue condition. To understand it better, we first note that the right hand 
side is exponentially small on account of the square of the Gamow factor e~ r °. Thus, 
ignoring for the moment the possibility that cos A c may vanish, either both e' and e" must 
be of order e~ v ° , or at the other extreme, one must be of order unity, and the other of 
order e~ 2r ° . Suppose that e' = 0(e~ 2r °), and e" = 0(1). Let us take A' = 1. Then from 
Eqs. ( |2T33|) and ( ggg ), we see that Q' = 0(1), while R' = 0(e~ 2rG ). It then follows from 



Eqs. ( p. 55 ) and fl2.56j ) that Q" ~ R" = 0(e r<3 ), and in turn from the double primed analogs 



of Eqs. ( gg3| ) and (^3§ that A" = 0(e" rG ). Lastly, since T" < T G , Eq. implies that 

B" is also exponentially small |[21|| . Hence, it follows that the entire wavefunction on the 
right hand side of the well, is exponentially small compared to the left hand part C . 
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In other words, there is negligible mixing of the states in the left and right hand well. This 
is exactly what we expect when the energies of the two states in the absence of tunneling 
differ by much more than the tunneling matrix element itself. 

The more interesting case, therefore, is that in which both e' and e" are of order e~ Ta . 
In the defining equations for the a's and /3's, Eqs. ( |2.33| ) and ( p.34[ ), we can then neglect 
the e's to very good approximation. This yields 



a' 1 /„, i\£'+h 



(3' V2n£'\ 
2tt' 



(2.59) 



where g n is the standard curvature correction f22|,|23fl in the phase integral expression for the 
tunnel splitting (see e.g., Eq. (4.10) of paper I). This quantity tends to 1 rapidly as n gets 
large: g = (n/e) 1 / 2 « 1.075, g x « 1.028, g 2 « 1.017, . . .. The right hand side of Eq. (p38|) 
can thus be written as A 2 (£',£")/ '4, where 

A{£\£") = -(^^) 1/2 (^o-o;o + ) 1/2 e- rG cosA c . (2.60) 

7T 

We further define 

e = |(e' + e") = E - \ + E + + {£' + + (f + ±> 0+ ) , (2.61) 

5 = e" - e' = (E_ -E + + (£' + |)w - - (f + |)wo+) . (2.62) 

With these definitions, Eq. ( |2.58| ) can be rewritten as 

e = ±-[5 2 + A 2 (£', £")} 1/2 . (2.63) 
2 



Equation (|2.63|) is the complete formal solution to the problem of tunneling in an asym- 



metric double well. Along with Eqs. (|2.60|) - (|2.62|) , it is the analog of the general phase 



integral formula for the tunnel splitting in a symmetric double well that we found in I [see 
Eq. (4.38) there]. Since there is no great need for having the final answer in simple closed 
form for the specfic problem (|1.1|) , and since the general procedure is fully explained in 
Sec. IV. E of I, we do not bother to extract the singular In J parts of the Tq integral. 

It is immediately obvious that the eigenvalues in Eq. ( |2.63| ) are what we would get from 
a two level system Hamiltonian 

nThS -2[A(£',£") -5 J' (2 ' 64) 

which is of course, just what we would expect. The quantity e is the energy measured from 
a convenient reference point, while 5, which depends on the fields h x , h z , and the quantum 
numbers £' and £" of the states whose mixing is being examined, is the offset between these 
energy levels in the absence of tunneling. Equation Q2.60 ) gives the tunneling amplitude 



between these levels when the offset is small, i.e., when the two levels are in approximate 
degeneracy. Note that although this amplitude is defined even for relatively large offsets — 
offsets comparable to the intrawell spacings ujq± — and indeed is not very sensitive to the 
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value of the offset, the concept of tunneling is physically sensible and useful only when the 
offset is comparable to or less than the amplitude A. If 5 ^> A, we get e ~ ±5/2, i.e., 
e" m 5, e' « A 2 / 8, or the other way around. Then by the argument given after Eq. ( [2.58 ), 
the mixing between the wells is negligible. 

One may wonder if the above conclusion does not invalidate the entire calculation. After 
all, we defined the shifts e' and e" assuming the wells were parabolic. Surely the corrections 
to the energies from cubic and higher order corrections to the potential are far larger than 
A. Since the offset between the levels must be tuned with exponential sensitivity, should not 
we know the energies of the levels before tunneling to the same sensitivity? The answer is 
no. The reason can be seen from Eq. ( |2.27 ). As we have seen, the amplitude of the growing 



part of the wavefunction plays a key role in the tunneling. From Eq. (|2.27|) , we see that 
this amplitude is only linearly dependent on e'. Hence, a small error in locating the absolute 
position of the level has little effect on the computed value of A. To say it another way, even 
though we know that the bottoms of the wells must be tuned to exponential accuracy to 
get significant mixing between the two wells, we cannot and need not determine the center 
of gravity of the two levels, e, to the same accuracy to determine the tunneling amolitude 
itself. This feature is also present in the symmetric case studied in I, although there it is 
not so apparent, since Herring's formula gives A directly without making reference to the 
absolute energy level. 

Secondly, it should be noted that Eq. (|2.63|) is nothing but the equation for the diabolo. 



The splitting vanishes only when 5 and A both vanish, which furnish the two conditions 
required to determine the diabolical point. Since both 5 and A will in general have linear 
terms in the deviation from the diabolical point, the energy surface is a double cone as 
asserted earlier. 

Thirdly, let us ask if we recover the results of paper I in the symmetric case, i.e., when 
h z = 0. Then E + = £L, M + = M_, and luq- = ujo+. The only sensible case is £' = i" = n, so 
that e' = e", and 5 = 0. The splitting is (up to an irrelevant sign) A(n, n), which is precisely 
the tunnel splitting A n computed in paper I. In addition, however, we now have more 
explicit information about the wavefunction. Proceeding as before, we see that Q' ~ A' ~ 1, 
R' ~ e~ r<3 , and B' ~ e _ ( rG_r '). (The double primed quantities are equal to their single 
primed counterparts.) The conclusion about B' is totally consistent with the approach 
in paper I, which was based on Herring's formula. There one takes the wavefunctions as 
(C m> d ± C- mt d) I v2, where G m ^ is the wavefunction of a state localized in the left well, and 
which decays away from that well in both directions. If we equate C m ^ with C' m la (and, 
therefore, C- m ^ with la ) in the central region, then we see that B', which by Eq. ( |2.16|) 
gives the magnitude of the growing part of C' m at m = m' c , is also the order of magnitude of 
C^ la , the decaying part of at m = m' c . 

I. What is the origin of the diabolical points? 

Lastly, it is extremely instructive to examine the problem when A = 0, i.e., cosA c = 0, 
without necessarily imposing the condition 5 = 0, for this gives insight into what causes the 
quenching of the tunnel splitting. Taking the imaginary parts of Eqs. ( |2.51| ) and ( |2.52| ) we 
see directly that we must have R' = R" = 0, and that 
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Q" = ±B'e r ° +r ", Q' = ±B"e r ° +r '. 



(2.65) 



Going back to Eq. (gjg) , we see that R' = requires either e' = or A' = 0, and likewise 
for R", e", and A". If 5 ^ 0, then both e' and e" can not vanish, and the only solution 
is e' 7^ 0, e" = 0, A" = Q" = B' = 0, (or the one obtained by interchanging single and 
double primes.) The only non-zero coefficients are A', Q', and B" . From Eqs. ( |2.33| ) and 
( p.65| ), we see that only one of these coefficients is independent, which can only be fixed by 
normalization. Thus, we see that the part C'^ 2 proportional to B" should really be regarded 
as the extreme right-hand tail of a state localized in the left well. If we denote this state by 
\L), and the wavefunction (m\L) by C^(m), then [15], 



C L {m) 



A'\vi(m)\ x l 2 exp [i J m qi(m')dm') , m < m' a , 

A'(— l) e D^(—z), m m' a tomR! m' b , 

<5 / |fi(m)| _1//2 exp i f™, qi{m')dm'^j , m' b < m < m b , (2.66) 

B"\vi{m)\' 1 / 2 exp (-i f™„ qi(m')dm'^j , m' b ' < m. 



[Note that in the second line, we wrote D^, not D u /, and that in the third line, we have 
not bothered to write the oscillatory exponential continuation in the region m' c < m < m!' c 
correctly — see Eq. ( ]2.39| ) — since our aim now is merely to indicate the general structure.] 
We can define a right-hand function Cn(m) analogously. Indeed it is now clear that the two 
state Hamiltonian ( |2.64| ) is a truncation of the full Hamiltonian (|1 . 1| ) in the \L), \R) basis. 

The above argument shows that for any H z (5^0), we can tune H x so that A vanishes, 
at which point, the energy eigenfunction is like Ci(m) [or Cn(m)], which is localized in 
one well and does not "see" the other well at all! In ordinary one dimensional quantum 
mechanics, this is of course impossible, since a wavefunction like Ci{m) which continues 
decaying with increasing m in the classically allowed region of the right hand well has the 
wrong sign of the curvature in that well. In fact, this argument does not depend on having 
5^0. If 5 = in addition to A = 0, C^m) and C^{m) are independent solutions of 
Schrodingers equation, as is any linear combination, since they are degenerate. 

The above point of view helps elucidate the origin of the quenching more clearly. In- 
deed, it is better to think about the non-symmetric situation (H z ^ 0) than the symmetric 
one (H z = 0). In continuum problems with a symmetric double well, Herring's formula 
gives the splitting as proportional to \ijjd{x){dip' d l 'dx)] x= o, where ipd(x) is a left-well localized 
wavefunction p^] . It is tempting to think that the splitting in the spin problem vanishes 
because the oscillatory envelope in C m) d in the central region allows the discrete analog of 
ipd{x) or ip'd( x ) to vanish at the midpoint. This reasoning is false, as one can see from a close 
examination of the symmetric case wavefunction, or even more clearly, by looking at the 
situation when A = but 5^0. The condition A = can not be reduced to a local property 
of the wavefunction such as its value or its slope at a particular point. Rather it is the global 
property that the phase integral A c be an odd integer times 7r/2. From this perspective, the 
quenching is perhaps better visualized as a manifestation of interfering Feynman trajectories 
and the Berry phase, even though the value of this phase is more easily found using the DPI 
method. 
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III. APPLICATION TO Fe 



Let us now apply our general formalism to Fe§. The problem of greatest interest is the 
location of the diabolical points, and for that we need only solve the conditions 5 = A = 0. 
We have already given formulas for the matrix elements and the critical curves in Eqs. (|2.5|) - 
( P.10| ). The problem that remains is to use these formulas to find S and A(£',£"). To keep 
the problem tractable and obtain answers in closed form, we will assume that h z is small. 
Specifically, we will assume that the reduced field h z defined above Eq. ( |2.8| ) is formally of 
order 1/ J. This enables us to evaluate the turning points and action integrals as expansions 
in powers of h z . Also, it is convenient to carry out all calculations in terms of the reduced 
variable \x. 

The first step is to obtain 5. For this, we need to analyze the critical curve Uo(fi). Its 
minima /iq± are found to be located at 



fi ± = ±yUo + 



+ 0{h% 



where 



to = (1 - hi) 112 . 

The quantities E±, u ±, and M± defined through Eq. (|2.18|) are given by 



E± = £/oGu ± 



J 



1 



2Xhl 




+ o(hl 



(3.1) 
(3.2) 

(3.3) 
(3.4) 
(3.5) 



Substituting these results in Eq. fl2.62|) , we obtain 



d(h z ,£',£") = 4v h z + 



where 



J 



d(h a 



V\h z 
J 



(£' + £" + l) Cl (h x ) + 0(J- s ) 



l-hl + \(l + 2hl) 
A(l - hi) 



(3.6) 



(3.7) 



The next step is to evaluate A(£',£"), or more precisely A c , the imaginary part of the 
phase integral defined in Eq. ( |2.50|) , since that by itself locates the diabolical points. This in 
turn requires expressions for x(m) and the points m' c and m". To obtain x( m )> we return to 
the Hamilton- Jacobi equation Eq. ( |2.4j ), write q = n + ix, and separate the equation into its 
real and imaginary parts. Eliminating K{m) from the two equations that result, we obtain 
a single equation for x( m )> which can be written as 
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At 2 (m)X 2 - g{m)X + 



t\{m) 
4t 2 (m) 



0; 



g(m) 



w(m) 



2U.(m) — E, 



(3.8) 
(3.9) 



where X = cosx(^)- [We can find the equation obeyed by K(m) similarly, and we discover 
that it is again Eq. ( |3.8| ) with X = cosh/t(m).] 

What value of E should we use in Eq. (|3T8| )? As stated earlier, the value of A is relatively 
insensitive to small changes in the absolute position of E. Thus we certainly needn't incor- 
porate the exponentially small shifts caused by the tunneling itself. Secondly, the tunneling 
is relevant only when 5 ~ A. Thus it suffices to set both e' and e" to in this part of the 
calculation. Then, Eqs. ( |2.61| ), ( |3.3| ), and (|3.4| ), we get 



E 



E_ + E + + {£' + | Vo_ + (f + |)cu 0H 

^VOfsi , on i i\ , ml J-2\ 



J 



+ £" + !) + 0(r 



Substituting this along with the formulas for w(m), ti(m), and t2(rn) in Eqs. 
yields the following equation for X: 



(3.10) 
and (p|), 



;i - A)(i - /i 2 )x 2 - [i + hi -f- c(ri\x + 



h 2 



1 - A 



0; 



CO*) 



^(£" -£') + l , (£' + £"+!)}. 



(3.11) 
(3.12) 



We have separated out the term in Eq. ( |3.11| ) as it is of order J -1 relative to the other 
terms. We can solve the quadratic equation and expand the result as a power series in (. 
Recalling that X = cos 2 \i we obtain 



cosx(yu) 



h. 



(1-A)(1-^) 



1 + 



CM 



2(1 - hi 



+ 0(C 2 



(3.13) 



By setting cosx = 1, we can determine the points [i' c and fi", for which it is now more 
convenient to write /i c ± instead. If the correction ( were absent, it is easy to see that these 
points would be at ±/i c0 , where 

= [(1 - A - h 2 x )/(l - A)] 1 / 2 . (3.14) 

Now, with ( ^ 0, we can find \i c ± as a series in J -1 . Up to leading corrections, we get 

1 



fi c± = ±/i 



cO 



2\/A/i c0 J 
The phase integral A c is given by 



A C = J 



x((J>)dfJ>- 



(3.15) 



(3.16) 



l-'c- 
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We will only evaluate this accurate to terms of order J°. Let us first consider the corrections 
entailed by replacing the limits by ±/U c0 . It is not difficult to show that without the ( term 
in Eq. (CT ), 



I 'I 2\ 1 /2 



(3.17) 



Inclusion of the £ correction does not change the square root approach to zero, and can not 
change the coefficient to leading order in J . From Eq. fl3.15|) , fi c ± differs from ±/x c o by 
terms of order J" 1 , so ignoring these shifts in the limits of the integral causes an error of 
order J -3 / 2 in the integral. Hence, we have 



I'cO 



Ac- J I x(/i)^ + 0(J- 1/2 ). 



(3.18) 



The error is smaller than 0(J°), and so may be ignored. The remaining integral may be 
done exactly as in paper I. We write the solution to Eq. ( [3.131 ) by x — Xo + A%, where xo 
is the solution when the £(/z) correction is ignored, and A% is the 0(() correction. The xo 
term can be integrated by parts, and yields 



A 



cO 



2J 

7TJ[1 



A) 



-1/2] 



(3.19) 



The Ax term yields 



AA r . = -J 



"McO 



C(a*) 



22 1 



hi 



V 2 ) 



dfj, 



-(£' + £" + I) 



71 



(3.20) 



where the result follows after an elementary integration by substitution. Note that the first 
term in £(/i) integrates to as it is odd. The total result for A c is (reverting to unsealed 
variables) 



A, 



A c0 + AA C 

7T 



2 J - (I' + £") - 2 J 



(3.21) 



It is now elementary to find the conditions for a diabolical point. The offset 5 becomes 
zero at a certain value of h z , and the tunneling amplitude A becomes zero when A c is an 
odd multiple of tv/2. Adding arguments £' and t" to indicate the level numbers becoming 
degenerate, the diabolicity conditions are, 



H z {£\£") V\{£"-£') 



2J 



1 + 



H x (£',£") y/T^X 



J 



J 



n 



4/^0 J 



(£' + £" + l) Cl (h x )+0(J- 2 ) 



1) 



n 



0,1, 



2 J - {£' 



1). 



(3.22) 
(3.23) 
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Note that what appears in these equations is J, not J. Secondly, the restrictions on the 
integer n in Eq. (|3.23|) are obtained, as explained in I, by demanding that A c be positive. 

We obtain the perfect centered rectangular lattice of diabolical points introduced in Sec. | 
if we ignore the C\ term in Eq. ( [3.22j ), and also the restriction H z /H c <C 1 used to perform 
the calculation. From the viewpoint of our DPI calculations, this exactness is somewhat 
mysterious. However, by the duality argument of Sec. [J if formulas ( p.22|) and (|3.22 ) are 
correct for small H z and large H x , then they also yield diabolical points for large H z and 
small H x . Of course if the former set of points corresponds to low lying levels, i.e., small 
values of £' and i", to which our analysis applies, the latter set of points corresponds to rather 
highly excited levels, to which the analysis does not apply prima facie. It is nevertheless 
surprising that the formulas should fit together so neatly. It is also somewhat surprising 
that the experimentally determined Fes diabolical points [for the cases (£', £") = (0, 0), (0, 1), 
and (0, 2)] should lie on a centered rectangular structure so closely, since, as is known, the 
value of the measured H x period on the H z = line is almost 50% different from that 
predicted by Eq. Q3.23| ). This feature is at present understood only on the basis of numerical 
diagonalization of the spin Hamiltonian including fourth order terms in J. An analytic 
approach to this aspect of the problem remains for the future. 
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FIGURES 

FIG. 1. Diabolical points for the Feg model Hamiltonian ( |1.1| ) for (a) J = 7/2, and (b) J = 4, 
as per the perfect lattice hypothesis. The numbers next to the light dashed lines show the numbers 
of levels that are simultaneously degenerate at all diabolical points on that line. By scaling the 
axes and reflecting about the heavy dashed line, we obtain points that are dual to one another; 
points on the line are self-dual. The pattern in the other three quadrants is obtained by reflecting 
about the H x and H z axes. 

FIG. 2. Critical energy curves for the Hamiltonian ([1 . 1|) . 

FIG. 3. Asymmetric double for a massive particle in one dimension, showing the turning points 
for an energy E. 

FIG. 4. Analog of Fig. ^ for the discrete case, showing only the left hand well. The critical 
curve U n is also not shown, as it does not determine any turning points for low lying levels. 
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TABLES 



TABLE I. Numerical test of perfect lattice hypothesis for J = 5. All diabolical points are 
obtainable from those listed here by duality and the symmetries H x — ► —H x , H z — ► —H z . In 
columns 4 and 5, i x = 2( J — n) — (£' + + 1), and i z = t" — £'. The next two columns give, for each 
value of A, the differences 8 X = \i x — 2JH x /H c ^Jl — A|, and S z = \i z — 2JH z /H c y/\\, where H x and 
H z are the numerically computed coordinates of the diabolical point. The numbers in parentheses 
give the power of 10 multiplying the number preceding. Where the error is given as 0.00, it is less 
than our machine accuracy. 
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